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DIFFUSION 
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The Wright-Fisher family of diffusion processes is a widely used 
class of evolutionary models. However, simulation is difficult because 
there is no known closed-form formula for its transition function. In 
this article we demonstrate that it is in fact possible to simulate 
exactly from a broad class of Wright-Fisher diffusion processes and 
their bridges. For those diffusions corresponding to reversible, neutral 
evolution, our key idea is to exploit an eigenfunction expansion of 
the transition function; this approach even applies to its infinite¬ 
dimensional analogue, the Fleming-Viot process. We then develop 
an exact rejection algorithm for processes with more general drift 
functions, including those modelling natural selection, using ideas 
from retrospective simulation. Our approach also yields methods for 
exact simulation of the moment dual of the Wright-Fisher diffusion, 
the ancestral process of an infinite-leaf Kingman coalescent tree. We 
believe our new perspective on diffusion simulation holds promise for 
other models admitting a transition eigenfunction expansion. 


1. Introduction. Monte Carlo simulation of diffusion processes is of 
great interest, as it underlies methods of statistical inference from discrete 
observations in a variety of applications (e.g. Golightly and Wilkinson, 2006, 
2008; Chib, Pitt and Shephard, 2010; Bladt and Sprensen, 2014; Bladt, Finch 
and Sprensen, 2016). Our interest in this paper is in the Wright-Fisher diffu¬ 
sion. This process is widely used for inference, especially in genetics, where 
it serves as a model for the evolution of the frequency Xt G [0,1] of a ge¬ 
netic variant, or allele, in a large randomly mating population. If there are 
two alternative alleles then the diffusion obeys a one-dimensional stochastic 
differential equation (sde) 

(1) dXt = -f{Xt)dt + VW(1 - Xt)dBt, Xo = xo,tG [0, T]. 
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The drift coefficient, 7 : [0,1] —)• M, can encompass a variety of evoiutionary 
forces. For example, 7 = /? where 

(2) /3(x) = ^[^ 1(1 - x) - 92 x\ + crx(l - x)[x + h{l - 2x)], 

describes a process with recurrent mutation between the two alleles, gov¬ 
erned by parameters 9i and 02 , and with (diploid) natural selection causing 
fitness differences between individuals with different numbers of copies of the 
allele, governed by parameters a and h. There is much interest among geneti¬ 
cists in inference from this and related diffusions (e.g. Williamson et ah, 2005; 
Bollback, York and Nielsen, 2008; Gutenkunst et ah, 2009; Malaspinas et ah, 
2012), and in the characteristics of the trajectories themselves (Schraiber, 
Griffiths and Evans, 2013; Zhao et ah, 2013), as discretely observed data 
are becoming more and more available (for example, as genetic time series 
data for ancient human DNA and for viral evolution within hosts). Beyond 
genetics, Wright-Fisher diffusions have been considered for applications in 
several other fields. For example, in finance they have been used as models 
for time-evolving regime probability, discount coefficients or state price (e.g. 
Delbaen and Shirakawa, 2002; Gourieroux and Jasiak, 2006); they have been 
proposed in biophysics as a model for ion channel dynamics (Dangerfield, 
Kay and Burrage, 2010; Dangerfield et ah, 2012); they have been studied as 
hidden Markov signals in filtering problems (Chaleyat-Maurel and Genon- 
Catalot, 2009; Papaspiliopoulos and Ruggiero, 2014; Papaspiliopoulos, Rug¬ 
giero and Spano, 2014); and in Bayesian statistics they have been proposed 
as models for time-evolving priors (Walker, Hatjispyros and Nicoleris, 2007; 
Favaro, Ruggiero and Walker, 2009; Griffiths and Spano, 2013; Mena and 
Ruggiero, 2016). 

Simulation from equation (1) is highly nontrivial because there is no 
known closed form expression for the transition function of the diffusion, 
even in the simple case 7 (x) = 0. In the absence of a method of exact 
simulation, it is necessary to turn to approximate alternatives such as an 
Euler-Maruyama scheme. Standard Euler-type methods fail here because 
simulated paths can easily leave the state space [ 0 , 1 ], and moreover stan¬ 
dard assumptions for weak and strong convergence typically require that 
the diffusion coefficient be Lipschitz continuous (see Kloeden and Platen, 
1999). Consequently, a number of specialized time-discretization methods 
have been developed for the Wright-Fisher diffusion with various drifts; when 
the drift is of the form of /3(x), see Schurz (1996) for 0i = 02 = o' = 0, Dan¬ 
gerfield et al. (2012) for a = 0, 0 i ,02 > 0, Schraiber, Griffiths and Evans 
(2013) for 01 = 02 = 0, /i = 1/2; and Neuenkirch and Szpruch (2014) 
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for cj = 0, 01,02 > 1- Other approaches include truncating a spectral ex¬ 
pansion of the transition function (Lukic, Hey and Chen, 2011; Song and 
Steinriicken, 2012; Steinriicken, Wang and Song, 2013), or numerical solu¬ 
tions of the Kolmogorov equations (Williamson et ah, 2005; Bollback, York 
and Nielsen, 2008; Gutenkunst et ah, 2009). The error introduced by these 
methods can be difficult to quantify and must often be tested empirically. 

In this paper we develop a novel and exact method for simulating the 
Wright-Fisher diffusion with general drift, as well as the corresponding 
bridges. By “exact” we mean that samples from the hnite-dimensional dis¬ 
tributions of the target diffusion can be recovered (up to the precision of a 
computer) without any approximation error. We build up our algorithm in 
stages, addressing how to simulate exactly from each of the following: 

1. The neutral Wright-Fisher diffusion; that is, with drift 

(3) a(x) = ^[0i(l - x) - 02x], 

where 0i,02 > 0 (Section 2). 

2. The neutral Wright-Fisher bridge (Section 3); informally, this is the 
process 

{^t)te[o,T] \ = y- 

3. The Wright-Fisher diffusion and its bridges with a very general class 
of drift functions (defined later), including drift /3(x) of the form (2) 
when 01,02 > 0 (Section 5). 

To achieve step 1 the key approach is to exploit an eigenfunction expansion 
representation of the transition function (see Griffiths and Spano, 2010, for 
review). The expansion admits a probabilistic interpretation and therefore 
lends itself to simulation techniques, but these techniques are not straight¬ 
forward to implement because the distributions involved are known only in 
infinite series form. Despite this hurdle, here we show that is possible to 
perform such simulation without error. The technique is very general and so 
we develop this section not just for the Wright-Fisher diffusion but for the 
Fleming- Viot process, its infinite-dimensional generalization. 

To achieve step 2 we obtain a new probabilistic description of the tran¬ 
sition function of a neutral Wright-Fisher bridge. This is again complicated 
by the appearance of distributions known only in infinite series form, but 
from which (we show) realizations can still be obtained by evaluating only 
a finite number of terms in the series. 

Finally, we generalize these techniques to nonneutral processes in step 
3. The eigenfunction expansion for a nonneutral process (Barbour, Ethier 
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and Griffiths, 2000) is probabilistically intractable, so we take a different 
approach: we use the simulated neutral processes as candidates in a rejec¬ 
tion algorithm. This uses a retrospective approach similar to that of the 
“exact algorithms” of Beskos and Roberts (2005), Beskos, Papaspiliopoulos 
and Roberts (2006, 2008), and Pollock, Johansen and Roberts (2016), which 
can return exact samples from a class of diffusions using Brownian motion 
as the candidate for rejection. We defer a detailed description to Section 5.1, 
but for now we note that a direct application of these techniques would re¬ 
quire that the target diffusion satisfy a number of regularity conditions, the 
most stringent perhaps that its law be absolutely continuous with respect 
to Brownian motion. The Wright-Fisher diffusion (1) fails in this regard, 
first because of its nonunit diffusion coefficient and second because of its 
finite boundaries. Although the first problem is easily solved via a Lamperti 
transformation [also known as Fisher’s transformation when applied to (1)], 
it is not clear how to deal with the second. Beskos, Papaspiliopoulos and 
Roberts (2008) point out that their exact algorithm can be adopted to the 
case of two finite entrance boundaries, but this approach requires a further 
technical condition [(A3) below] which does not always hold here. When it 
does hold, this approach becomes arbitrarily inefficient when the diffusion 
is proximate to the boundary (Jenkins, 2013); in any case, the boundaries 
of the Wright-Fisher diffusion can be exit, regular reflecting, or entrance, 
depending on the parameters of f3(x). But now we are armed with the abil¬ 
ity to simulate the neutral Wright-Fisher process, which serves as a far 
more promising candidate than Brownian motion in a rejection algorithm; 
specifically, it is known that the law of a nonneutral process is absolutely 
continuous with respect to its neutral counterpart (Dawson, 1978; Ethier 
and Kurtz, 1993). We develop these ideas in full in Section 5. We also re¬ 
mark that a related approach is taken by Schraiber, Griffiths and Evans 
(2013), who were interested in simulating nonneutral Wright-Fisher bridges 
in the absence of mutation. In this context one can condition each sample 
path to remain in (0,1) (otherwise the path could be absorbed at 0 or 1 and 
could not terminate at a pre-specified point), rendering the boundaries in¬ 
accessible. They show that the appropriate candidate in this case is a Bessel 
process of dimension 4, whose boundary at 0 is also of entrance type. How¬ 
ever, their method is not exact in the sense given above, since the rejection 
probabilities are approximated via numerical integration. Furthermore, the 
Radon-Nikodym derivative of the Wright-Fisher process with respect to the 
Bessel(4) process is not bounded (another approach developed in Jenkins 
(2013) for a single entrance boundary suffers a similar limitation). In any 
case a direct comparison with our method is not possible since here we tackle 
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01 , 6*2 > 0 rather than 6 i = 62 = 0 . 

The remainder of the paper is structured as follows. In Section 4 we 
discuss practical considerations of the algorithm; Section 6 fills in one last 
gap by showing how to construct a nonneutral Wright-Fisher bridge; Section 
7 discusses extensions of the algorithm; and Section 8 contains the proofs of 
our results. 

2. Simulating the neutral Wright-Fisher process. In this section 
we demonstrate how exact simulation from the neutral Wright-Fisher dif¬ 
fusion can be achieved. To aid the exposition we first focus on a one¬ 
dimensional process, and then later extend this idea to the Fleming-Viot 
process. 

2 . 1 . A transition density expansion in one dimension. Consider the dif¬ 
fusion satisfying (1) with drift (3). Denote its law by and its transition 

density f{x, •; t)- Throughout this paper we assume 0i, 02 > 0; then f{x, ■; t) 
is a probability density. We will exploit the following probabilistic represen¬ 
tation of the transition function’s eigenfunction expansion (Griffiths, 1979; 
Tavare, 1984; Ethier and Griffiths, 1993; Griffiths and Spano, 2010): 

00 m 

(4) fix, y,t) = Y^ q^{t) I 3 m^^il)Ve^+i,e 2 +m-iiy), 

m=0 1=0 

where 

Bm,xil) = (^^x\l - / = 0,1,... ,m, 

is the probability mass function (pmf) of a binomial random variable, 

is the probability density function (pdf) of a beta random variable, 0 = 
01 -|- 02 , and {qfnit) ■ m = 0,1 ,...} are the transition functions of a certain 
death process A^it) with an entrance boundary at 00 . More precisely, let 
{A^{t) : t > 0} be a pure death process on N such that ^^(0) = n almost 
surely and whose only transitions are m 1 —?• m — 1 at rate m{m + 6 — 1 )/ 2 , 
for each m = 1, 2,..., n. Then q^it) = limn^oo IP’(^n(^) = ''^)- 

The representation (4) has a natural interpretation in terms of Kingman’s 
coalescent, which is the moment dual to the Wright-Fisher diffusion. The 
ancestral process A^dt) represents the number of lineages surviving a time 
t back in an infinite-leaf coalescent tree, when lineages are lost both by coa¬ 
lescence and by mutation. For our purposes, the mixture expression (4) also 
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Algorithm 1: Simulating from the transition density /(x, •;t) of the 
neutral Wright-Fisher diffusion with mutation. 

1 Simulate 

2 Given A^{t) = m, simulate L ~ Binomial(m, x). 

3 Given L — I, simulate Y ~ Beta(@i + 1,02 + m — 1). 

4 return Y. 


provides an immediate method for simulating from f{x, ■; t). We summarize 
this idea in Algorithm 1, which first appeared in Griffiths and Li (1983). 

Steps 2 and 3 of Algorithm 1 are straightforward, but Step 1 requires the 
PMF {qm{t) : m = 0,1,...}, which is not available in closed form. Griffiths 
and Li (1983) used a numerical approximation, but in the following subsec¬ 
tion we show it is in fact possible to simulate from this distribution without 
error. 


2.2. Simulating the ancestral process of Kingman’s coalescent. Our goal 
in this subsection is to obtain exact samples from the discrete random vari¬ 
able with PMF {q^it) : m = 0,1,...}. Were each qf^it) available in closed 
form, then standard inversion sampling would return exact samples from this 
distribution (see for example Devroye, 1986, Gh. 2): for U ~ Uniform[0,1], 

{ M 

M G N : ql{t) > U 

m=0 

is distributed according to {q^it) ■ m = 0,1,...}. However, q^it) is known 
only as an infinite series (Griffiths, 1980; Tavare, 1984): 


Qm{t)= ""“Le where 

(5) k=m 

0 _ {0 + ‘i.k — 1){9 + m)(fc_i) 

“ m\{k-m)\ ■ 

Here we have used the notation 0(3.) := r(a-|-x)/r(a) for a > 0 and x > — 1. 

Despite the apparently infinite amount of computation needed to evaluate 
(5), we now show that it is nonetheless possible to return exact samples from 
this distribution by a variant of the alternating series method (Devroye, 
1986, Ch. 4), which we summarize for a discrete random variable A on N as 
follows. Suppose X has pmf {pm : m = 0,1,...} of the form 


00 

Pm = 

k=0 
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and such that 

( 6 ) bk{m) 4 , 0 as A: —)• 00 , for each m. 

Then for each M, K 

M 2K+1 M M 2K 

T-{M) =Y.Y. < Ep- ^ E 

m=0 k=0 m=0 m=0 k=0 

Furthermore, T^{M) t P(X < M) and T^{M) 4 , P(X < M) as iF —)■ 00 . 
Hence, for U ~ Uniform[0,1] and 

Ko{M) := inf {K £N:T-{M)>U or r+(M) < U} , 

the quantity inf |m G N : | can be computed from finitely 

many terms and is exactly distributed according to {pm '■ m = 0,1,...}. 

This approach can be applied—with some modification—to {g^(t) : m = 
0 , 1 ,...} by the following proposition, which says that the required condition 

( 6 ) holds with the possible exception of the first few terms in m. For those 
exceptional terms, ( 6 ) still holds beyond the hrst few terms in k, and there 
is an easy way to check when this condition has been reached. 

Proposition 1. Let , the relevant coeffi¬ 

cient in (5), and let 

(7) := inf {i > 0 : . 

Then 

(i) Cm^'^ < 00 , for all m. 

(a) b^^’^\m) 4 , 0 as k ^ 00 for all k >m + Cm^\ and 
(in) = 0 for all m > Dq’^\ where for e G [ 0 , 1 ), 

( 8 ) := inf jz > Q - V 0 : (0 + 2 Z + l)e-™ < 1 - e| . 

(The parameter e is introduced for later use.) As a result of Proposition 1, 
we need only to make the following adjustment to the alternating series 
method: If m < then precompute terms in qf^{t) until the hrst time 

that the coefficients in (5) begin to decay; we know by Proposition l(ii) that 
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Algorithm 2: Simulating from the ancestral process Aoo(t) of King- 
man’s coalescent with mutation. _ 

1 Set m <— 0, ko <— 0, k <— (fco). 

2 Simulate U ~ Uniform[0,1]- 

3 repeat 

4 I Set km i — [C^’®V2l [eq. (7)]. 

5 while 5^ (m) < U < (m) do 

6 I Set k i —fc + (1,1,..., 1) 

7 end 

8 

9 if {m) > U then 

10 I return m 

11 else if S'^(m) < U then 

12 Set k <— (fco, fci, • • ■, fcm,0). 

13 Set m — m + 1. 

14 end 

15 until false 


this decay then continues indefinitely. To allow for the number of computed 
coefficients to depend on m we introduce k = {ko, ki,, kM) and 

(9) 

M 2km+t M 2km 

S+{M) := Y.'- 

m=0 i=0 m=0 i=0 

We summarize this procedure in Algorithm 2. 

Of course, the false condition in line 15 of Algorithm 2 is never met, but 
Proposition 1 guarantees that the algorithm still halts in finite time. Further 
performance considerations of this algorithm are discussed in Section 4. Also 
note that computed coefficients in (5) can also be stored for future calls 
to this algorithm. 

2.3. A transition density expansion in higher dimensions. It is worth 
pointing out that an interesting by-product of Proposition 1 (and of Algo¬ 
rithm 2) is the possibility of simulating exactly from the transition function 
of the (parent-independent) neutral Wright-Fisher diffusion in any dimen¬ 
sion, even in infinite dimensions. Wright-Fisher diffusions in d dimensions 
can be seen as d-dimensional projections of a so-called neutral (parent- 
independent) Fleming-Viot measure-valued diffusion = {fj,t ■ t > 0) with 
state space Aii{E), the set of all the probability measures on a given (Polish) 
type space E, equipped with the Borel sigma-algebra induced by the weak 
convergence topology. Given a total mutation parameter 6 and a mutation 
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distribution Pq £ A4i(E), the process fi is reversible with stationary distri¬ 
bution given by the Dirichlet process with parameter {6 ,Pq), here denoted 
with n^ pp, characterized by Dirichlet finite-dimensional distributions: 


n, 


e,Po 




) G dxi, } oc 


\2 = 1 




_ 2=1 


(^ 1 ? ■ ■ • ? ^d) 


for any d and every measurable partition Ai,..., ^4^ of where = 

{(xi,..., Xrf) E [0, Xi = 1}. 

The transition function describing the evolution of ^ admits a probabilis¬ 
tic series expansion as mixture of (posterior) Dirichlet processes: 


( 10 ) 


oo « 

p{fl,du-t) = / p^^{dCl,...,d^m)nQ ep{du), 

^ ^ I pm ’ 0+m ' 9 + m 

m=0 

t > 0,/i, zz E 


where denotes the n-fold /r-product measure and r]m ■= YllLi % 
(see Ethier and Griffiths, 1993). The coefficients of the series expansion are 
given by i.i.d. samples (the ^-random variables) from the starting measure, 
/i, of random size given by the coalescent lines-of-descent counting process 
with distribution An algorithm for simulating from the tran¬ 

sition function (10) is thus the following modification of Algorithm 1. 


Algorithm 3: Simulating from the transition density p(/r, ■ ;t) of the 
neutral Fleming-Viot process with parent-independent mutation. 

1 Simulate A^{t). 

2 Given A^{t) = m, simulate Ci; ■ ■ • > Cm 

3 Given = dm, simulate ~ ^ e 

4 return v. 


Notice that step 3 requires sampling a (potentially infinite-dimensional) 
random measure distributed according to a Dirichlet process. Techniques for 
exact sampling from a Dirichlet process have been available in the literature 
[e.g. Papaspiliopoulos and Roberts (2008) and Walker (2007)] for quite some 
time. Hence Algorithm 2 provides a way of filling the only missing gap 
(Step 1 of Algorithm 3) to make the transition function (10) viable for exact 
simulation. When E consists of d points {d E N) the process reduces to the 
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(d — l)-dimensional Wright-Fisher diffusion, thus Algorithm 1 is viable for 
exact simulation of neutral {d — l)-dimensional extensions of the Wright- 
Fisher diffusion (1) with drift (3). 

3. Simulating a neutral Wright-Fisher bridge. In this section we 
demonstrate how exact simulation from the neutral Wright-Fisher diffusion 
bridges can be achieved, via a new probabilistic description of its transition 
density. For the remainder of the paper we return to processes of dimension 
one. 


3.1. A transition density expansion. Section 2 provides a method for 
returning exact samples from f{x,-]t) for any fixed x G [0,1] and t > 0. 
The density of a point y G (0,1) at time s in a Wright-Fisher bridge from 
X at time 0 to 2; at time t is given by (Fitzsimmons, Pitman and Yor, 1992; 
Schraiber, Griffiths and Evans, 2013): 


( 11 ) 


fz,t{x,y]s) 


f{x,y;s)f{y,z;t- s) 
f{x,z-,t) 


0 < s < t, 


with /(•,•;■) as in (4). Motivated by (4), our goal is to facilitate easy simu¬ 
lation from fz,tix,y;s) by putting it into mixture form. For the rest of this 
section we assume 0 < x,y, z < 1. 


Proposition 2. The transition density of a Wright-Fisher bridge has 
expansion 

oc 00 m k 

(12) f^^t{x,y;s) = EEEEiTr 

m=0 k=0 1=0 j=0 


where 

(13) 


Pm,k,l,j — ^m,x\l)F>ei+j,e2+k-j{z)VAAQ.^^j^i^Q.^j^rn-l\k\3) f(^xZ't) ' 
for 0 < I < m and 0 < j < k, and P^’ki j’^^ ~ 1* otherwise, where 


F>Ma,b-,k{j) = 


k\ B{a + j,b + k - j) 
3 ) B{a,b) 


is the PMF of a beta-binomial random variable on {0,1,..., /c}. 
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By Proposition 2, we recognize equation (12) as a mixture of beta-distributed 
random variables, with mixture weights defining a PMF ^ ^ 

N} on Thus, the following algorithm returns exact samples from fz,t{x, y, s). 


Algorithm 4: Simulating from the transition density fz^t{x,y;s) of a 
bridge of the neutral Wright-Fisher diffusion with mutation. 

1 Simulate {M,K,L,J) ~ ^ N- (13)]. 

2 Given (M, K, L, J) = simulate 

Y ~ Beta(0i + l + j,92 + m + k — I — j). 

3 return Y 


Again, while step 2 of Algorithm 4 is straightforward. Step 1 is compli¬ 
cated by the appearance of ~ s)/f{x,z]t) in (13); each term in 

this ratio is known only as an infinite series, as we have seen. We address 
this in the following subsection. 

3.2. Simulating from the discrete random variable on with pmf 
iPm’kl’j’^'^ • ^ a-PPly the alternating series approach 

of Section 2.2 separately to each of qfnis), q'^{t — s), and f{x, z; t), and then 
combine these to obtain monotonically converging upper and lower bounds 
on (13). The first two terms have been dealt with already in Section 2.2, so 
it remains to take a similar approach for f{x, z\ t). Note that this problem— 
the pointwise evaluation of f{x,z;t )—is separate from (and in this case, 
harder than) actually simulating from f{x, • ;^)- 

To employ the alternating series approach, use (4) and (5) to write 

OO OO 

(14) 

m=0 k=m 

where 

and Lm ~ Binomial(m, x). Our strategy is to group the triangular array of 
coefficients (4*^’*’^^) k>m in such a way that, with the exception of the first 
few terms, they exhibit a property analogous to (6). We will compare terms 
in the sequence (di)i=o,i,... ol antidiagonals, defined by 

m 

/ . nt — u, 1, . . . , 

i=0 i=0 


(15) 


J _ 

U2m / ^ ^m+j,m—j ’ 


^2m+l 
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(see Figure 2), and dropping the superscript for convenience. Notice that the 
coefficients within each entry of this sequence are all multiplied by the same 
sign in (14), so that f{x, z-,t) = (Iq —di + d 2 — d 3 + ... will be our alternating 
sequence. The main complication in this approach is to find explicitly the 
first i for which the coefficients (d*) begin to decrease in magnitude. To this 
end we define 

(16) ;= inf |m > 0 : 2j > for all j = 0 ,..., m| , 

which simply provides the first entry in {d 2 m) for which every member of 
the corresponding antidiagonal is decaying as a function of its first index. 
We now need the following lemma. 

Lemma 1. Let Lm ~ Binomial{m, x) and 

(17) := ^ + 

z 1 — z 

Then for all m G N, 

(18) E['D 0 ^+Lrr,+l,e 2 +m+l-Lm + liz)] < , 62 +^ 1 -L^i^)] ■ 

Using Lemma 1 we are in a position to obtain the required analogue of 
property (6). 

Proposition 3. Let (di)j=o,i,..., , and be defined as 

in (8), (15), (16), and (17), respectively, and e G (0,1). Then 

d2m+2 ^ d2m+l < d2m 

for all m > V V /e. 

We can now combine Propositions 1 and 3 in order to construct a sequence 
amenable to simulation from ^ (13)], via the 

alternating series method. 


Proposition 4. 


p(s^t,9,x,z) 

and 



Define 

:= V V V V /e, 




Lj=o 


m] 


iAt-sfi) 


J Lj=o 


(fc) 






i=0 


X 
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Then for 


+ < em,fe,;j(2y + 3) < < era,uj(2v + 2) < em,kui^P)‘ 


In other words, for sufficiently large v the odd and even terms in the 
sequence {em,k,i,j{v))^=Q provide monotonically converging lower and upper 
bounds on respectively, and “sufficiently large” can be verified 

explicitly. 

The above results are summarized in Algorithm 5. To explore we 
introduce for convenience a bijective pairing function S : N —?• N^, such that 
S(n) = {m,k,l,j). As in Algorithm 2, we also introduce v = {vo,vi, ... ,vn) 
and 

N N 

K7(^) := + '^), ■= eS(n)(2Un). 

n=0 n=Q 


Algorithm 5: Simulating from the discrete random variable on with 
pmf : rn,k,l,j E N}. 

1 Set n <— 0, Vo <— 0, v i — (vq). 

2 Simulate U ~ Uniform[0,1]. 

3 repeat 

Set Vn i — [Ts(„)/2] [eq. (19)]. 
while V~{n) <U < V^{n) do 
Set V i — u + (1,1,..., 1) 

end 


7 

8 
9 

10 

11 

12 

13 

14 


if V~{n) > U then 
return S(n) 
else if V^{n) < U then 

Set V <— (i;o,ui,... ,Un,0). 
Set n i — n + 1. 
end 


15 until false 


4. Performance of algorithms for neutral processes. There are 
several easy improvements to the underlying Algorithm 2. For example, we 
are free to vary the order of inspection of each m by any finite permutation 
of N, and we found a dramatic improvement by radiating outwards from (an 
approximation of) the mode of '■ m = 0,1 ,...} than to start at m = 0 
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and work upwards. Our approximation used in Theorem 1 below. It 

may also be possible to improve on Algorithm 2 by allowing different q^it) 
to be refined at different rates, i.e. by using a vector other than (1,1 ,..., 1) 
in Step 6; we do not explore that here. 

A crucial quantity governing the efficiency of our algorithms is the num- 
ber of coefficients 6^’ (m) we must compute in Algorithm 2; these in turn 
depend on the quantities and Cm^'^ (recall Proposition 1). These quan¬ 
tities are in general manageably small, suggesting that the number of coeffi¬ 
cients that need to be computed in Algorithms 2 and 5 (line 4 in each) should 
not be excessive. One exception to this observation is when t is very small, 
for which the number of relevant coefficients grows quickly. The following 
result makes precise the complexity of Algorithm 2. 

Proposition 5. As t ^ 0, 

(i) = O(t-i). 

(ii) rnaxmCm^^ = 0{t~^ log{t~^)). 

(Hi) for any k > 0 . 

Let denote the total number of eoefficients that must be computed in 

an implementation of Algorithm 2. Then < oo, and in particular 

(iv) E[A'(hO] = o{t~^^^^^), for any k > 0 . 

The growth in Algorithm 2 as t 0 is closely related to the well known 
numerical instability of (5) as t —)• 0 (Griffiths, 1984), which afflicts any 
method based on the expansion (4) (or an expansion using a basis of orthog¬ 
onal polynomials, which is equivalent to (5); Griffiths and Spano, 2010). In 
any practical implementation of our algorithm we are obliged to use an ap¬ 
proximation should the separation between two points be very small, and we 
found Algorithm 2 to fail for t < 0.05 or so. One option is to revert to the 
Euler-Maruyama approximation for small t. Alternatively, there has been 
much previous work in coalescent theory on approximating the distribution 
(5) (e.g. Griffiths, 1984, 2006; Jewett and Rosenberg, 2014); by inserting 
those approximations into Algorithm 2 they readily define new algorithms 
for approximate simulation of the dual diffusion. We consider the following 
result, due to Griffiths (1984, Theorem 4). 
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Theorem 1 (Griffiths (1984)). Suppose (3 = ^{9 — l)t, and let 
= (^ M ))2 


V 

T/ien P ((j4^(t) —< x) —^ 4)(x) as t —)• 0, where <!)(•) is 
the cumulative distribution function (CDpj of a standard normal random 
variable. 

[The statement in Griffiths (1984) is missing the factor To apply this 
approximation in practice when t is small, we replace line 1 in Algorithm 1 
with 

1’ Simulate Aoo(t) ~ N{p!d’^\ (cr^*’®))^) and round it to the nearest non¬ 
negative integer. 

4.1. Comparison with Euler-Maruyama simulation. To check the cor¬ 
rectness of our algorithm and to compare its performance to Euler-Maruyama 
simulation, we performed the following experiment. We fixed 6 i = 62 = 1/2, 
explored various fixed values of x and t, and for each parameter combi¬ 
nation drew 10,000 samples from /(x,-;t) using Algorithm 1. To quan¬ 
tify whether the resulting sample was consistent with the true distribution 
f{x,-;t), we performed a one-sample Kolmogorov-Smirnov (K-S) test. We 
then performed the same experiment instead using Euler simulation with 
various stepsizes 6 . For this purpose we used the Balanced Implicit Split 
Step (BISS) algorithm of Dangerfield et al. (2012), an Euler-type algorithm 
with some advanced modifications that guarantee each sample path stays 
within [0,1], and which is state-of-the-art for 61,82 > 0. (Their algorithm 
has an additional tuning parameter e; we followed their recommendation 
and set e = 6/4.) 

To obtain an accurate expression for the CDF of the true distribution 
for use within the K-S statistic, we exploited the fact that, in the special 
case 61 = 02 = 1/2, a Lamperti transformation of (1) (or conversion to 
Stratonovich form) leads to 

Xt = ^(1 - cos Bt), 

where {Bt)t>o is a Brownian motion commenced from arccos(l — 2x) and 
reflecting at 0 and vr. A rapidly converging series expression for the CDF of 


^(?? +/?)^ (^1 + — 2?7^/3 /3 7^ 0, 



where 


/3 = 0 , 


/ 3 / 0 , 
13 = 0. 
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Bt (and hence of Xt) is available (Linetsky, 2005, eq. (26)); the hrst 1000 
terms in the series suggested convergence to machine precision and were 
used as the reference cdf. 

As is evident from Table 1, exact simulation strongly outperforms the 
BISS algorithm over almost all the parameter combinations investigated. 
Over a timescale of t ^ 0.1, errors in the Euler-type method accumulate suf¬ 
ficiently fast that the resulting samples are easy to reject in a K-S test. Even 
reducing the stepsize so that its running time is several orders of magnitude 
greater than the exact method provides only a modest improvement to the 
quality of the sample. Note also that the performance of the BISS algorithm 
deteriorates for paths started close to the boundary (compare x = 0.01 with 
X = 0.5), whereas the exact method is indifferent to starting position. One 
reason p-values are persistently small for the BISS algorithm is that sample 
paths are constrained by construction to stay inside [<5/4,1 — (5/4], yet the 
narrow region close to the boundaries is precisely where we expect to hnd 
much of the probability mass for many choices of parameters. An example 
of how this affects the resulting transition density is given in Eigure 1. By 
contrast, in no application of the K-S test to samples from our method would 
we have rejected at level 0.05 the hypothesis that they were drawn from the 
true distribution. Only over short timescales and away from the boundaries 
(Table 1, t < 0.05 and x = 0.5) do the two methods seem comparable. At 
t = 0.05, the same computational investment as in the exact method but 
applied to the BISS method buys a stepsize of about 6 = 10“^^, which is ade¬ 
quate in the interior of the state space (x = 0.5) but not near the boundaries 
{x = 0.01). 

When t > 0.05 it is also possible to compare the exact method both 
with and without the approximation of Griffiths (1984) (see pl5); the two 
versions exhibit similar running times and generate bona fide samples from 
the true distribution (according to a K-S test) for moderate t. However, the 
approximate version deteriorates (as reported by the K-S p-value) for large 
t (see Table 1, t = 5), away from its asymptotic regime. Thus a suitable 
rule-of-thumb is to use the exact algorithm for t > 0.05 and its approximate 
version for t < 0.05, with the two methods in good agreement in their region 
of overlap in t. 

We investigated several other choices of x and t with predictable results 
(not shown); for example, performance of the exact method seems unrelated 
to starting position x, while the BISS method deteriorates even further as 

X —> 0. 
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Table 1 

Comparison of exact simulation methods for the neutral Wright-Fisher diffusion. BISS; 

Algorithm of Dangerfield et al. (2012) with stepsize S. Exact; Algorithm 1. Exact’; 
Algorithm 1 with the approximation of Griffiths (1984) described on pl5. Tabulated are 
the computing time needed to simulate 10,000 sample paths and the p-value of a K-S test 
applied to the resulting collection of endpoints. Paths are initiated at Xq = x and run for 
length t. Mutation parameters are 9i = 62 = 1/2. 





t = 

0.01 




X 

= 0.01 



X = 0.5 


Method 

<5 

Time (s) 

p-value 

Method 

S Time (s) 

p-value 


10 -^ 

0.03 

< 10-36U 


10-5 002 

8.0 X 10-6 


10-"* 

0.05 

< 10-36° 


10-"* 0.05 

0.36 

BISS 

10-® 

0.30 

< 10-3°° 

BISS 

10-5 0.30 

0.22 


10-® 

2.83 

< 10-3°° 


10-6 2.78 

0.30 


10-^ 

27.61 

< 10-3°° 


10-'3 27.44 

0.78 

Exact’ 

- 

0.19 

0.94 

Exact’ 

0.17 

0.18 




t = 

0.05 




X 

= 0.01 



X = 0.5 


Method 

<5 

Time (s) 

p-value 

Method 

8 Time (s) 

p-value 


10 -^ 

0.04 

< 10-366 


10-5 004 

9.0 X 10-3 


10-"* 

0.16 

< 10-3°° 


10-"* 0.16 

0.35 

BISS 

10“® 

1.40 

< 10-3°° 

BISS 

10-5 1.39 

0.53 


10-® 

13.77 

< 10-3°° 


10-6 13.66 

0.02 


10-^ 

138.08 

< 10-3°° 


10-'3 137.06 

0.72 

Exact 

- 

0.35 

0.09 

Exact 

0.34 

0.64 

Exact’ 

- 

0.17 

0.30 

Exact’ 

0.17 

0.97 




t = 

: 0.5 




X 

= 0.01 



X = 0.5 


Method 

<5 

Time (s) 

p-value 

Method 

5 Time (s) 

p-value 


10“^ 

0.16 

< 10 - 3 °° 


10-6 al6 

1.2 X lO-^’* 


IQ-'* 

1.43 

< 10 - 3 °° 


10-3 1.43 

9.0 X 10-38 

BISS 

10-^5 

14.07 

< 10 - 3 °° 

BISS 

10-5 14.08 

9.4 X 10-33' 


10-® 

137.82 

< 10 - 3 °° 


10-6 138.23 

6.3 X 10-33 


10-^ 

1378.33 

< 10 - 3 °° 


lO-’’ 1368.33 

8.1 X 10-33 

Exact 

- 

0.19 

0.16 

Exact 

0.19 

0.81 

Exact’ 

- 

0.17 

0.21 

Exact’ 

0.17 

0.60 




t - 

= 5 




X 

= 0.01 



X = 0.5 


Method 

s 

Time (s) 

p-value 

Method 

5 Time (s) 

p-value 



0.17 

< 10 - 3 ™ 


10 -^ ore 

< IO-36'J 


10-3 

1.43 

< 10-3°0 


10-6 1.43 

< 10-3°° 

BISS 

lO-'^ 

14.01 

< 10-3°° 

BISS 

10-"3 14.04 

< 10-3°° 


10-5 

137.95 

< 10-3°° 


10-5 138.09 

4.7 X 10-30° 


10-6 

1375.75 

< 10-3°° 


10-6 1378.54 

1.6 X 10-30° 

Exact 

- 

0.18 

0.58 

Exact 

0.18 

0.88 

Exact’ 

- 

0.17 

1.2 X 10-'‘3’ 

Exact’ 

0.17 

0.49 
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B, 


(a) 

Fig 1. Histogram of 10,000 simulated points using (a) Algorithm 1 and (b) the Euler- 
type BISS algorithm of Dangerfield et al. (2012). Parameters are 61=62 = 1/2, x = 0.5, 
t = 0.5, 5 = 10“®. For visual clarity, plotted are samples of the driving reflecting Brownian 
motion Bt = arccos(l — 2Xt) rather than Xt (since the density function of Bt—shown as 
a solid line—is bounded and can be calculated; see main text). 



5. Simulating the nonneutral Wright-Fisher process. In this sec¬ 
tion we develop an exact rejection algorithm for simulating from the Wright- 
Fisher diffusion (1) with general drift. We make use of retrospective sampling 
techniques for the exact simulation of diffusions, which we first summarize. 

5.1. Overview of the exact algorithm. Here we give a brief overview of the 
exact algorithm (EA) of Beskos and Roberts (2005), Beskos, Papaspiliopou- 
los and Roberts (2006, 2008), and Pollock, Johansen and Roberts (2016), 
and we refer the reader to those papers for further details. The EA returns 
a recipe for simulating the sample paths of a diffusion X = {Xt)t^[o^T] sat¬ 
isfying the SDE 

(20) dXt = fl{Xt)dt + dBt, Xq = xo,t € [0,T], 

with fl assumed to satisfy the requirements for (20) to admit a unique, weak 
solution. Denote the law of such a process, our target, by Qa:o- The idea of 
the EA is to use Brownian motion started at xq, whose law will be denoted 
'WxQ, as the candidate process in a rejection sampling algorithm. The goal 
is then to write down the rejection probability, which is possible under the 
following assumptions: 

(Al) The Radon-Nikodym derivative of with respect to Wa,,, exists and 
is given by Girsanov’s formula, 

^ (X) = exp {£ ,(W - 1 . 


(21) 
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(A 2 ) ^ G 

(A3) (j){x) := is bounded below, by </>“ say. 

(A4) A{x) := fj,{z)dz is bounded above, by A+ say. 

Using (A1-A4) and Ito’s lemma, (21) can be re-expressed as 

( 22 ) ^^(X) oc exp{A(Xr) - A+}exp|-y [4>{Xt) - (j)~]dt 

Written in this form, the right hand side of (22) is less than or equal 
to 1, and therefore provides the required rejection probability. To make 
the accept/reject decision, it is necessary to construct an event occurring 
with probability (22). This is easy to achieve given a realized sample path 
(Xt)jg[Qj-] ~ Wa;^, but obtaining such a path would require an infinite 
amount of computation. Instead, note that the right hand term in (22) is 
the probability that all points in a Poisson point process $ = {{tj,ipj) : j = 
0,1,...} of unit rate on [0, T] x [0, oo) lie in the epigraph of 11 —)• [(p{Xt) — 4>~], 
and this event can be determined by simulating X only at times ti,t 2 , ■ ■ ■■ 
Thus, the following algorithm returns a (random) collection of skeleton 
points from X 



Algorithm 6: Exact algorithm for simulating the path of a diffusion 
process with law 


1 repeat 

2 I Simulate a Poisson point process on [0,T] x [0,oo). 

3 Simulate U ~ Uniform[0,1]. 

4 Given $ = {{tj,ipj) : j = 0,1,...}, simulate B ~ Wa, at times 
(tj)j=o,i,... and at time T. 

5 if — 4>~ < ipj for all j = 0,1,... and U < exp{A{BT) — A"*"} 

then 


return {{tj,BtA : j = 0,1 ...} U {(T, 5^)}. 


end 


8 until false 


Once a skeleton has been accepted, further points can be filled in as re¬ 
quired by sampling from Brownian bridges; no further reference to Qa, is 
necessary. If (p is bounded above, by <p'^ say, then Algorithm 6 can be imple¬ 
mented with finite computation by thinning $ to a Poisson Point process 
on [0,T] X [0,4>^ — 4>~]', |$| is then almost surely finite. [However, this re¬ 
quirement on cj) can be relaxed (Beskos, Papaspiliopoulos and Roberts, 2006, 
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2008).] 

We remark that assumption (A4) is restrictive, and can in fact be removed 
by using a certain biased Brownian motion as an alternative candidate; this 
also improves the efficiency of the algorithm (Beskos and Roberts, 2005). We 
present the EA in the form above since an analogue of (A4) does hold for the 
Wright-Fisher diffusion, and in any case a ‘biased Wright-Fisher diffusion’ 
is not available. 


5.2. Exact algorithm for the Wright-Fisher diffusion. As noted earlier, 
the requirements (A1“A4) need not hold when our target is the Wright- 
Fisher diffusion. However, related techniques can be used if the candidate 
process is chosen to be another Wright-Fisher diffusion, and in particular 
one with the same mutation parameters. Denote the target law WF.y_3;(j [with 
drift 7], and the candidate law WFQ_a:o [with drift (3)]. We will write 7 G WW 
if there exists a drift a. of the form (3) such that the following hold: 

(WFl) The Radon-Nikodym derivative of with respect to WFa,a;o ex¬ 

ists and is given by Girsanov’s formula. 


(23) 


dWF. 


7 : 2:0 


dWF, 


(A) = exp 


a,a:o 


7(W) - a{Xt) 


dXt 


Xtil - Xt) 

2 Jo Xt{l - Xt) 


dt 


(WF2) 7 is continuously differentiable on (0,1). 

(WF3) (p{x) is bounded on (0,1): (l)~ < 4>{x) < where 


fix) 




x(l 


a^{x) 

x) 


+ 7'(a;) 


a'{x) - [7(x) 


a{x)] 


1 - 2 x 
x(l — x) 


(WF4) A{x) := dz is bounded above, by A+ say. 

Specific conditions on a and 7 to satisfy (WFl) are detailed e.g. in Theorem 
8.6.8 in 0ksendal (2003). (This theorem imposes some unduly restrictive 
conditions to ensure that the SDE has a unique weak solution, but this can 
be established for (1) by other means; Pardoux, 2009, Ch. 4.) Following 
Section 5.1, we apply Ito’s lemma to A{x) to re-express (23) as 


We recognize the rightmost term in (24) as the probability that no points 
in a Poisson point process on [0, T] x [0, — (/>“] lie in the epigraph of 
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Algorithm 7: Exact algorithm for simulating the path of a diffusion 
process with law 

1 repeat 

2 Simulate a Poisson point process on [0,r] x 

3 Simulate U ~ Uniform [0,1]. 

4 Given '■ j = 0, U • • • , J}, simulate X ~ WFc^a; at times 

and at time T. 

5 if — 4>~ < ipi for all j = 0,1,..., J and U < exp{A(AT) — A^} then 

6 I return : j = 0,1 •••,>/} U {(T, At)}. 

7 end 

8 until false 


1 1 —)■ 4>{Xt)—4>~. Hence, Algorithm 7 returns exact samples from Step 

4 of Algorithm 7 is achieved via Algorithm 1. Once a skeleton is accepted, 
further points can be filled in via Algorithm 4. 

5.3. A class of drifts for which exact simulation is possible. The class 
ysfT defines the set of drifts for which exact simulation is possible. Here 
we show that WT contains all the most popular population genetics diffu¬ 
sion processes with mutation and selection (including frequency-dependent 
selection), whose drift admits the general form 

(25) ''y{x) = a{x) + x{l — x)ri{x), 

where a is as in (3) for some 0 i ,02 > 0 and ry is a reasonably regular 
function of x. The case of diploid selection (2) is recovered by setting rj(x) oc 
(x -|- h{l — 2x)). For general r] the properties of diffusion processes with 
drift (25) and of the corresponding genealogies are studied, among others, 
by Coop and Griffiths (2004). 

Proposition 6. Any drift of the form (25), for a as in (3) for some 
01,6*2 > 0 and for rj continuously differentiable in (0,1), is a member ofVdT. 

It is worth noting that the additive structure in the drift 7 (that is, with 
a selection component added to the mutation component a) is a widely ac¬ 
cepted and theoretically justified property of all population genetics diffusion 
models (see e.g. the discussion in Karlin and Taylor, 1981, pl86-187). 

5.3.1. Example: Wright-Fisher diffusion with diploid selection. By Propo¬ 
sition 6 , the drift fd [eq. (2)] satisfies fd G YdfF . In fact, in the haploid case 
{h = 1 / 2 ), there is much simplification; f is a quadratic polynomial for 
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Table 2 

Performance of Algorithm 1 applied to the Wright-Fisher diffusion with symmetric 
mutation and additive selection. Each row reports means (per accepted path) across a 
simulation to generate 1,000 accepted paths. Paths were initiated at Xq = x and run for 
time t, with mutation parameters 6 i — 62 = 0.01 and selection parameter a (and 
h — 0.5^. Reported are the total numbers (per accepted path) of: attempts, Poisson points 
simulated, coefficients needed, random variables generated (that is, the aggregate 

of all draws from underlying constituent distributions: uniforms, betas, and so on), the 
number of times the simulation resorted to the approximation of Theorem 1, and the 

running time in milliseconds. 


(7=1 


t 

X 

Attempts 

Poisson 

points 

Coefficients 

Random 

variables 

G1984 

Time (ms) 

0.1 

0.5 

1.27 

0.00 

116.07 

7.37 

0.00 

0.007 

0.1 

0.25 

1.40 

0.00 

132.87 

8.01 

0.00 

0.008 

0.1 

0.01 

1.58 

0.01 

141.34 

8.91 

0.01 

0.009 

0.5 

0.5 

1.23 

0.03 

18.65 

7.22 

0.00 

0.003 

0.5 

0.25 

1.48 

0.03 

21.60 

8.45 

0.01 

0.003 

0.5 

0.01 

1.58 

0.03 

23.84 

9.02 

0.01 

0.003 

5.0 

0.5 

1.29 

0.24 

5.23 

8.58 

0.00 

0.002 

5.0 

0.25 

1.49 

0.27 

6.18 

9.64 

0.01 

0.002 

5.0 

0.01 

1.67 

0.28 

7.36 

10.79 

0.01 

0.002 




a 

= 10 




0.1 

0.5 

11.83 

3.72 

995.92 

62.68 

1.83 

0.062 

0.1 

0.25 

41.69 

12.97 

3714.82 

224.95 

7.86 

0.225 

0.1 

0.01 

145.73 

45.96 

14937.52 

856.21 

45.88 

0.879 

0.5 

0.5 

13.16 

20.96 

641.52 

109.00 

2.47 

0.054 

0.5 

0.25 

43.82 

69.05 

2729.80 

399.95 

10.89 

0.214 

0.5 

0.01 

149.21 

235.34 

17044.43 

1869.54 

71.54 

1.185 


which analytic bounds are available, and A(^x) = axj^. We implemented 
our exact algorithm on this model, and investigated its performance by con¬ 
sidering several combinations of parameters; results are shown in Table 2. 
For moderate selection (|(t| = 1) the algorithm is extremely efficient, with 
only slightly more than one candidate needed per acceptance. Furthermore, 
most simulations resulted in zero Poisson points. These results are quite ro¬ 
bust to the length of the path t, the initial frequency x, and the sign of the 
selection parameter. For stronger selection (a = 10) we observe some dete¬ 
rioration in efficiency because of the greater mismatch between candidate 
and target paths—to the extent that simulation of paths of length t = 5.0 
became prohibitive. Nonetheless, it is still feasible to simulate a collection 
of shorter paths in a few seconds (and to string these together to construct 
longer paths, if necessary). 
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To make this observation more precise, Beskos, Papaspiliopoulos and 
Roberts (2006, Proposition 3) obtained an explicit upper bound on the ex¬ 
pected number of Poisson points required of Algorithm 6, and hence on the 
computational complexity of the algorithm. A careful reading of their re¬ 
sult shows that it relies on the existence of the bounds cj)^ but does not 
depend on the laws of the diffusions involved, and carries over easily to the 
Wright-Fisher case. We therefore omit a proof of the following. 

Proposition 7. Let denote the number of Poisson points required 

until the first aceepted path. Then 

< (^+ - 

An immediate consequence of Proposition 7 is that the complexity of 
simulating a path of length KT is 0{{(j)'^ — )kt+a+'^ 

KT —> oo. However, superior performance can be achieved by splitting 
the problem into K separate simulations; the complexity is then — 

(j)~)KTe^'^^~‘^ )T+A+^^ which is linear in K as in Beskos, Papaspiliopoulos 
and Roberts (2006). 

As this is a statement about the complexity of the algorithm as the path 
length increases, it continues to hold even if we account for the increasing 
cost associated with each Poisson point as T —0, as quantified by Proposi¬ 
tion 5. In practice one might wish to optimize the choice of K and T for a 
given path length KT. In our application, we must be prepared to introduce 
an additional constraint in order to fix T some distance away from 0 (and, 
as we recommend above, the choice T > 0.05 seems adequate). 

6. Simulating a nonneutral Wright-Fisher bridge. For complete¬ 
ness, here we provide an algorithm for simulating a nonneutral Wright-Fisher 
bridge (Algorithm 8). This follows immediately from the previous sections; 
the only modification is to note that the appropriate candidate process is the 
corresponding neutral bridge, which can be simulated via Algorithm 4. The 
rest follows exactly as in the Brownian case; see Beskos, Papaspiliopoulos 
and Roberts (2006, Section 6.2) for details. 

7. Discussion. In this paper we have shown how to simulate exactly 
from the scalar Wright-Fisher diffusion, as well as a number of important 
and closely related processes: these include the ancestral process of an oo-leaf 
Kingman coalescent tree, the Fleming-Viot process with parent-independent 
mutation, the nonneutral Wright-Fisher diffusion, and neutral and nonneu¬ 
tral Wright-Fisher bridges. Some interesting open problems remain, includ¬ 
ing mutation operators which do not lead to reversible diffusions, and the 
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Algorithm 8: Exact algorithm for simulating the path of a diffusion 
process with law conditioned on Xt = y- 

1 repeat 

2 Simulate a Poisson point process on [0,r] x 

3 Given ■ j = 0, !> • • • i J}, simulate X ~ (WFc.x | Xt = y) at times 

4 if <l){Xt -) — < V'j for all j = 0,1, ..., J then 

5 I return {(G.XtJ : j = 0,1 U {(r,i/)}. 

6 end 

7 until false 


problem of sampling from (d — l)-dimensional (2 < d < oo) Wright-Fisher 
bridges. 

It is also important to notice that, in order to employ the machinery pro¬ 
posed in this paper, the mutation parameters 61,62 in the a-component of 
a general drift 7 need to be both positive: models of the form (25) with 
no mutation (a = 0 ) or with one-directional mutation (with only one mu¬ 
tation parameter positive and the other null) have at least one absorbing 
boundary, and therefore there cannot be absolute continuity with respect to 
a stationary Wright-Fisher diffusion with a transition density expansion of 
the form of (12). For such cases, series expansions are in fact available with 
a structure similar to (4) (see e.g. Ethier and Griffiths, 1993) and we believe 
it should be relatively simple to adapt our method to encompass selection 
models with absorbing boundaries, a goal we do not pursue here. For choices 
of drifts 7 more general than (25), arising possibly in models beyond pop¬ 
ulation genetics, it is harder to specify conditions for (WFl) verifiable in a 
simple way by inspection of 7 . The determination of which drift functions 7 
guarantee (assuming identical diffusion coefficients) absolute continuity with 
respect to a Wright-Fisher process, seems to be, to our knowledge, an open 
problem. We expect that most of what is affected by the drift pertains to the 
behaviour of the process at the boundaries. The problem, however, might 
be delicate and not just limited to making sure that 7 prevents the bound¬ 
aries from being absorbing: for example, if for the diffusion with drift 7 the 
boundaries are entrance or reflecting, the rate of escape from the bound¬ 
aries might still make its sample paths qualitatively different from those of 
any Wright-Fisher diffusion, even within the same boundary regime, respec¬ 
tively entrance or reflecting. To support this conjecture, we remark that the 
very same circumstances induce mutual singularity in squared Bessel pro¬ 
cesses (whose diffusion coefficient is ^/x hence quite similar near zero to the 
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Wright-Fisher diffusion): it is well-known that any two squared Bessel pro¬ 
cesses starting at 0 are mutually singular whenever their drifts differ, even if 
they are both within the same boundary regime (see Pitman and Yor, 1981, 
Lemma (2.1) and references therein). 

In a wider perspective, we believe that the approach proposed here might 
serve as a template for developing new techniques for sampling exactly from 
diffusion processes by means of non-Brownian bridges whose transition func¬ 
tion admits a transparent transition function expansion. 


8. Proofs. Proof of Proposition 1. First suppose m > 0. 

(i) Note that 


(26) 


um 

4+1 


dm 


(m) 

(m) 


6 + m + k — l 6 + 2k + l ( 2 k+e)t 

k — m + 1 6 + 2k — 1 


fm{k)e 


{2k-\-G)t 

2 


say. Treat fm{k) as having domain M; it then suffices to show that < 

0 for all sufficiently large k. Then the right hand side of (26) is subsequently 
decreasing in k monotonically to 0. Part (i) follows for the finite k {= m + i) 
for which the right hand side of (26) drops below 1. Routine calculations 
show that {fmYik) < 0 for all k > {^2{m — 1) + 6 — 0)/2. 

(ii) Note that 


^2{m-l)+e -e 

2 


< 



Ve-e 

2 


< 



1 

+ -<m, 


so in fact ifYiYi^) < 0 k > m, and thus as soon as {tti) < 

for some fc, it must also hold for all subsequent k. 

(iii) The right hand side of (26) can be made independent of m by noting 
that 


(27) ak)<fY{k) = e + 2k + l. 

(f 0^ 

Thus for Cm = 0 to hold we need m to exceed the upper of the two 
solutions on M of 

{9 + 2k + l)e"(2fc+0)42 ^ 

The definition of is one way to express this condition, since the max¬ 
imum of (6 + 2k which separates the two solutions is at 

k = [j — Finally, consider the special case m = 0. If 0 > 1 then sim¬ 

ilar arguments as in (i-ii) above continue to hold. However, if 0 < I then 
in fact (/o)^(fe) > 0 for all k, with /o(^) continuous on /c > 1. But then 
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/q(A:) < /o(c>o) = 1 for > 1, so /o(^)e (2^+^')V2 <; i fQj. all A: > 1 and 
hence (i-ii) still hold, with Cq’^^ <1. □ 

Proof of Proposition 2. This follows by substituting (4) into (11), multi¬ 
plying by B{6i + l+j,92 + m-l + k- j)/B{6i + l + j,92 + m - I + k - j), 
and rearranging. □ 


Proof of Lemma 1. First suppose I < \rnz\. Then, using r(x-|-l) = xr(x). 


]P(-bm+l = 0^01 

m + l , 9 + m , 


^{Lm — 0^ei+z,e2+m-;(^) 


< 


< 


m + 1 


m + 1 — mz 
1 — X 


(1-x)- 


9 + m 


1 - z 


92 + rn — mz 


(1-z) 




maximizing the term in square brackets first in I and then in m, noting for 
the last inequality that this term is increasing in m. Hence, summing over 
1 = 0 ,..., [mz], 


(28) 


[mz] 

^ ^ IP(Tm+l — /)Xl0j^-|_;^02-l-m+l—i (^) 
1=0 


< 


1 — X 


1-z 


\rnz\ 

Y^nB^ 


1=0 


l)Vg 

i+i,6»2+m-z(^)- 


By a similar argument, for I > [mz\: 

X 

^{Bm+1 = I + g2-|-m+l-(Z-|-l)(2:) < -F{Lm = O^0i-|-Z,6»2-l-m-Z(^), 

(this time it is crucial we compare T^+i = 1 + 1 with Lm = I, rather than 
with Lm = 1 + 1), and hence 


m+1 

(29) E 

1= [mz\ +1 

m 

< - Yh ^{Bm = O^0i+Z,02+m-K^)- 

l=\rnz\ 
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m 


Cm + l ,m + l 


C 





Cm+2,m 


CO.O 


Cl,l C2,l 


Cl,0 C2,0 


Cm +1 ,m — 1 C 7 Ti + 2 ,m — 1 1 




C2m,0-* -C2m + 1.0-<- C2m+2,0 


do di d2 da 


- r- 7 k 

d2m d2m + l d2m+2 


Fig 2. Computation of The sum of each antidiagonals (dashed) defines the se¬ 

quence (di)i=o,i,... ■ To show that di+i < di terms are paired off as shown by the arrows; 
that is, the coefficient at the head of a set of arrows is greater in magnitude than the sum 
of the terms at its tails. 


Finally, sum (28) and (29) to yield (18), noting that the overlap on the right- 
hand side at I = [mz\ necessitates the given definition of (instead of 

the simpler bound V |). □ 

Proof of Proposition 3. First, since m > we know j > for 

(f (f 0'\ 

j = 0,..., m and hence by Proposition 1 that — — j). 

Now multiply this inequality by ¥\Ve.^+Lrr,-jfi 2 +m-j-L.^_j{z)] to yield 


(30) 


(x,z,t,9) {x,z,t,e) 


Thus, summing over j = 0,1,..., m, 


(31) 


ix,z,t,9) 

/ J ^m-\-l-\-j,m—j 


3=0 


< 


E 


(x,z,t,9) 
^m+j,m-j > 


1=0 


which says precisely that < d 2 m- We also need to show that d 2 m -\-2 < 

d 2 m-\-i, but this case is more subtle since the left hand side is a sum over one 
more term than the right. Instead, we will increment the first index in (30) 







28 


P. A. JENKINS & D. SPANO 


and then sum over j = 1 ,... ,m: 


(32) 


E 


{x,z,t,e) 

^m+ 2 +j,m—j 




E {x,z,t, 0 ) 

^m+l+j,m-j • 


It then suffices to show 


(33) 


{x,z,t, 0 ) {x,z,t, 6 ) {x^z,t, 6 ) 

' ^771+2,m ’ 


for if we sum (32) and (33) we obtain d 2 m +2 < d 2 m+i as required (Figure 2). 
To show (33), first note that 


{x,z,t, 9 ) 

(x,z,t, 0 ) 

'"k,m 


4 *+iV) 


= fLik)e 


(2k + 9)t 

—^ < (0 + 2A: + l)e- 


(2k+9)t 

2 


with f^{k) defined as in (26) and the inequality following from (27). Hence, 
choosing k = m + 1 and noting that m > , 

/oA\ ^ I o; I (x,z,t, 0 ) . /, \ {x,z,t, 6 ) 

(34) cl^+ 2 ,J < {0 + 2 k + l)e 2 c+gj < (1 - 


Second, note that 

(x,z,t, 0 ) 

{x,z,t,e) 


(35) 


^ 1 6* + 2m E[V 0 ^+L^^^,e 2 +m+i-Lrr,+i (z)] 

m + ie + m E[V0^+Lm,e2+ni-L^iz)] 

^ 1 

- m + 1 HT^g^+L^,e2+m-LA^)] 

< e, 


using Lemma 1 and m + 1 > je for the last inequality. Rearrange 

(35) and sum with (34) to get (33). □ 


Proof of Proposition 4. This follows immediately from Propositions 1 and 
3. It can also be viewed as an application of Proposition 1 of Beskos, Pa- 
paspiliopoulos and Roberts (2008) to a function g{ui,U 2 ,U 3 ) oc U 1 U 2 /U 3 . □ 


Proof of Proposition 5. (i) Using (26) and (27), 


(36) 


dt,8) 

^K+1 


(m) 


(2K+e)t 


(2K+0)t 


dm 


(m) 


= f^{K)e -^ <(0 + 2K + l)e- 2 
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The right-hand side of (36) is less than 1 provided 

log(6' + 2m + l) 9 

K >-, 

t 2’ 

which implies that — ^ = 0{t~^) as t —)• 0. 

(ii) Continuing, 


(j 0 ) log( 6 * -L 2m + 1) ^0 log( 6 * -L 2D^ ’ ^ -F 1) _ 0 
^ t 2 t 2’ 

with the right-hand side independent of m and asymptotically 0{t~^ log(t“^)) 
as t —)• 0 by (hi) below. 

(f 0^ (f 

(hi) By inspection of the dehnition (8) of Dq ’ , in order to ensure K > Dq ’ 

\i suffices that 

IV log( 6 ' + 2K + 1) 

^ t ’ 

which holds for sufficiently small f if iV ~ for any hxed k > 0. 

Hence, = o(f“^^+'^^) as f —> 0. 

(iv) Parts (i-iii) cover those terms that must be precomputed in Algorithm 2: 
For the random number of remaining terms we may assume K > m + Cm , 
so that bfc+iim) < bK{m). Write Q%{t) := P(A^(f) < M). Our aim is to 
show that these random remaining terms do not add to the complexity of 
the calculation beyond (iii); we achieve this by determining the complexity 
of 


{ M oo 

K > max (m + ^ ^ (-l)^-™6fc(m) < (5 

^_Q 

as f —>• 0, for a fixed 5 > 0. This is the appropriate quantity to look at, 
since if iV > cf’^\M) then \Q%{t) — S^{M)\ < 6, when 2ki > K; i = 
0,..., M. Furthermore, averaging over the uniform random variable U drawn 
for inversion sampling, the total number of coefficients \ U = u needed 

to determine that < u < Qmi't) is given by the number of terms 

needed to bound both our estimates of Qm-i{t) arid Qm{t) away from u: 


E(Ard’^i) =E[E(iVi*’®i I U)] 


(37) 


< 


E 


rQlit) 


m=0 •^Qm-i(b 


^ Q?n(t) 


(m) + C 
- 1 /. '' ■' 


(t,0) 


(m — l)]dn. 


We will show that if iV ~ f (i+'^l as t —>• 0, for a fixed k > 0, then we can 
attain the stated growth condition on Using that the right-hand 
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(f 0) 

side of (36) is decreasing in K for K > m + Cm , for each () < 1 we can find 
a constant such that for K > Hence, 


(38) 


M oo M oo 

Y1 <YjY1 


m=0 k=K 


m=0k=K 


M oo 

< E E C"-'"&A(m) 

m=0 k=K 


M 


E 

m=0 


bKjrn) 

1-c ■ 


Routine calculations show that, for m > 1, 


e _ e {k - m){6 + m +k - 1) 
«fc,m+l - (m+l)(0 + m) 


^ e (k — 1)(^ + k) ^ n 

— ^km n/n I — ^kl 


{k-i){e + k) 

2(0 + 1 ) 


1 fc-i 


.(9) 


2(0 + 1 ) 

and thus, applying Stirling’s formula to oil ~ 

( 39 ) bK{m) = Kio^K-KiK+e-i)t/2 < 


for some constants C 2 ^\ (which again exist by our assumption about 

the asymptotic growth of K). In the special case m = 0, Stirling’s formula 

also yields a^g ~ /c+ and so the inequalities in (39) continue to hold. 
Combining (38) with (39) we find 


M oo 

m=0 k=K 


for some which is less than S if 
(40) A 


(This is not a tight bound but suffices in the calculations below.) In sum¬ 
mary, if K satisfies both K > and (40) then K > cf’^\M). Inte- 


log- i 

^cf(M + l)\ 

V t ^ 1 
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grating over 6: 


( 41 ) f 

Jc 


Ql{t) 








{m)du= f cf’^\m)d6 

Jo 




< 


= Q^it) 


^ + l/ 7 log 

tl + K 


c^f\m + 1 ) 


d6 


+ A/- 


log(c^(m + 1 )) + 1 - logg^(t) 


It remains to show that the resulting expression (41) can be summed over 
m, which is possible by Theorem 1; in particular, 


9m (i) = 


1 


y^27r((Th’®))^ 


exp 


(m — 
2 (cj(he ))2 


+ 0 ( 1 )) 


Hence, by Jensen’s inequality, 

00 

9m(^)log"i = IE[logH^(t)] < logE[H^(t)] = 0(logt“^), 


m=0 


Y 9m(i)[-log9m(i)] 


m=0 




1 . 


= O(logt-I) + -E[x 2 ] + 0 ( 1 ) = O(logt-i), 


where X ~ N{0, 1). Combining these results with (41), we obtain 


rQlit) 


Y /. = 0{t l^+''l) + 0(t ^/^logt ^) 


m=0 


= 0 (t-(i+^)), 


showing that the first term on the right of (37) is 0{t The sec¬ 

ond term is also 0(t“l^'’''’’l) by a similar argument. Since k was arbitrary, 
E[ivM)] = o(t-(i+'’’)). □ 

Proof of Proposition 6. For a diffusion with drift 7 , since 7 is continuous 
on ( 0 , 1 ), then 

£ = £ [d\Wt{l-Xt) + 2ijiXt)aiXt)] dt<oo, 
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then Novikov’s condition is satisfied and a Girsanov transform exists with 
respect to the neutral WF diffusion with drift a, i.e. (WFl) holds. In par¬ 
ticular, (23) reads 

expj^ r]{Xt) dXt - ^ [nHXt)Xt{l-Xt) + 2r,{Xt)aiXt)] dt 

The function r] is also continuously differentiable in (0,1), so ^'{x) — a'(x) = 
r]\x)x{l — x) -|- r]{x){l — 2x) is continuous and 

~ 1 

(43) 4>{x) = - [x(l — x){rf‘{x) + Vj' {x)) -|- 2?7(x)a(x)] 

which, being itself continuous in [0,1], is then bounded in (0,1), and (WF3) 
follows. 

(WF2) and (WF4) are obvious. Thus (42) has a version of the form (24) 
with (f) as in (43) and A{x) = ?](z) dz. □ 
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